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ABSTRACT 

The existence of supermassive black holes lurking in the centers of galaxies and of stellar binary 
systems containing a black hole with a few solar masses has been established beyond reasonable doubt. 
The idea that black holes of intermediate masses (~ 1000 M Q ) may exist in globular star clusters has 
gained credence over recent years but no conclusive evidence has been established yet. An attractive 
feature of this hypothesis is the potential to not only disrupt solar-type stars but also compact white 
dwarf stars. In close encounters the white dwarfs can be sufficiently compressed to thermonuclearly 
explode. The detection of an undcrluminous thermonuclear explosion accompanied by a soft, transient 
X-ray signal would be compelling evidence for the presence of intermediate mass black holes in stellar 
clusters. In this paper we focus on the numerical techniques used to simulate the entire disruption 
process from the initial parabolic orbit, over the nuclear energy release during tidal compression, the 
subsequent ejection of freshly synthesized material and the formation process of an accretion disk 
around the black hole. 

Subject headings: meshfree Lagrangian hydrodynamics, nuclear reactions, reactive flows, black holes 



1. INTRODUCTION 

The existence of two classes of black holes is well- 
established: supermassive black holes with masses be- 
yond 10 6 M Q are thought to lurk in the centers of most 
galaxies ( j23j) and black holes with just a few solar masses 
have been identified as the unseen components in X-ray 
binary systems (|12p. Black holes with Mbh ~ 1000 Mq, 
so-called "intermediate mass black holes" , represent a 
plausible, but so far still unconfirmed "missing link" be- 
tween these two well-established classes of black holes. 
In recent years, the stellar dynamics in the centers of 
some globular clusters has been interpreted as being the 
result the gravitational interaction with an intermedi- 
ate mass black hole(j5] [HI E3 0. Further circumstan- 
tial evidence comes from ultraluminous, compact X-ray 
sources in young star clusters Q35| I2ip and from n-body 
simulations (22) that indicate that runaway collisions in 
dense young star clusters can lead to rapidly growing 
black holes. None of these arguments is conclusive by 
itself, but their different nature suggests that the pos- 
sibility of intermediate mass black holes must be taken 
seriously. 

The tidal disruption of white dwarfs offers the unique 
possibility to explore the presence of an intermediate 
mass black hole. The corresponding disruption processes 
have been explored with various approximations in ear- 
lier studies (TTil [Ml [5J) , our simulations for the first time 
explore the full evolution from the initial parabolic orbit 
over the disruption process to the subsequent build-up 
of an accretion disk. We have performed a large set of 
calculations to identify observational signatures that can 
corroborate or, alternatively, rule out the existence of in- 
termediate mass black holes. In this paper we focus on 
the numerical techniques that have been employed in this 
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disruption study, a detailed discussion of the astrophys- 
ical implications will be given elsewhere (!27l |2"5| . 

2. NUMERICAL METHODS 

The simulation of a white dwarf disruption by a black 
hole needs to follow the gas dynamics from the initial 
spherical star through the distortion and compression 
while approaching the black hole to the subsequent ex- 
pansion phase and the formation of an accretion disk. 
Of paramount importance for the dynamical evolution is 
the inclusion of the feedback from the nuclear reactions 
that are triggered by the tidal compression. In the fol- 
lowing we will briefly sketch the methods employed in 
our simulations. 

2.1. Hydrodynamics 

Due to the highly variable geometry and the impor- 
tance of the strict numerical conservation of physically 
conserved quantities, we use the smoothed particle hy- 
drodynamics method (SPH) to discretize the equations 
of an ideal fluid. Using the SPH approximations (J3j ITS} 
the conservation of mass, momentum and energy trans- 
late into 
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Here p a is the density at the position of particle a, m b 
the (constant) mass of particle b, W ab — W(\r a — f%\, h ab ) 
the cubic spline kernel (|15|1 with compact support whose 
width is set by the average of the smoothing lengths of 
particle a and b, h ab — (h a + h b )/2. P a refers to the 

gas pressure, f a , g rav to the gravitational acceleration of 
particle a and Tl ab is the artificial viscosity tensor, see 
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below, and v ab = v a — Vf, with v being the particle ve- 
locity. The quantity u a is the specific internal energy of 
particle a, e n uc,a is the thermonuclear energy generation, 
calculated in a operator split fashion as described below. 
Since SPH is a Lagrangian method, and we ignore mix- 
ing, the compositional evolution is a local phenomenon. 
In cases where the geometry of the gas distribution varies 
substantially, it is advisable to adapt the local resolution, 
i.e. the smoothing length, according to the changes in the 
matter density. We achieve this by scaling the smoothing 
length with the density according to 
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where the index labels the quantities at the beginning of 
the simulation. The initial smoothing lengths are chosen 
so that each particle has 100 neighbors 4 . Taking the 
Lagrangian time derivative of both sides of Eq. Q and 
using the Lagrangian form of the continuity equation, 
one finds an evolution equation for the smoothing length 
of each particle: 
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This equation is integrated together with the other ODEs 
required for hydrodynamics, Eqs. (8j), and possi- 

ble changes in the abundances in the case ofnuclear reac- 
tions. During the integration the resulting new neighbor 
number is constantly monitored for each particle and, if 
necessary, an iteration of the smoothing length is per- 
formed to keep the neighbor number in the desired range 
of between 80 and 120. 

The SPH equations derived from a Lagrangian (|3"T1 [T7)l 
yield different symmetries in the particle indices and ad- 
ditional multiplicative factors with values close to unity, 
but a recent comparison (29) between these two sets of 
equations showed only very minor differences in practi- 
cal applications. 

We use the artificial viscosity tensor in the form given 
in ()16p . but with the following important modifications: 
i) the viscosity parameters a and (3 (commonly set to 
constants a = 1, f3 = 2) are replaced by a — > a ab = 
(a a + a b )(f a + /&)/4, (3 2a ab , where f k is the so- 
called Balsara-switch© to suppress spurious forces in 
pure shear flows, and ii) the viscosity parameter a k is 
determined by evolving an additional equation(19) with 
a decay and a source term(| 
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In the absence of shocks a decays to a m - m on a time 
scale T a = h a / {0.1 c a ), where c a is the sound velocity. In 
a shock a rises rapidly to capture the shock properly. A 

4 A "neighbor" is aparticle b that yields a non-zero contribution 
to sums in Eqs. |l}-|3|- 



further discussion of this form of artificial viscosity can 
be found in (j25jl . 

To demonstrate the ability of this scheme to capture 
shocks without spurious post-shock oscillations we show 
in Fig. [l]the results of a standard, one-dimensional Sod 
shock tube test(30) (with a polytropic gas of adiabatic 
exponent L = 1.4) at t=0.17. For this test we used 1000 
SPH particles, a max = 1.5 and a m i n = 0.1. The numeri- 
cal solution (circles) agrees excellently with the exact so- 
lution (solid line), the shock itself is spread over about 10 
particles. In the third panel ("pressure") we have over- 
laid the value of the viscosity parameter a. It deviates 
substantially from its minimum value only in the direct 
vicinity of the shock where it reaches values of about 1.3. 
This is to be compared to the "standard" values a = 1 
and (3=2 that are commonly applied in astrophysical 
simulations to each particle irrespective of the necessity 
to resolve a shock. This scheme has been shown (j2l)l l2"S"jl 
to minimize possible artifacts due to artificial viscosity. 
Note that no efforts have been made to include the effects 
of heat transport. Therefore, deflagration-type combus- 
tion cannot be handled adequately with the current code. 
As will be discussed below, deflagration is of no impor- 
tance for the disruption processes that we discuss here. 

2.2. Equation of state 

The system of fluid equations (HI , (p| and (|3]) needs to 
be closed by an equation of state (EOS) that is appropri- 
ate for white dwarf matter. We use the HELMHOLTZ 
EOS developed by the Center for Astrophysical Ther- 
monuclear Flashes at the University of Chicago. It ac- 
cepts an externally calculated nuclear composition which 
facilitates the coupling to reaction networks. The ions 
are treated as a Maxwell-Boltzmann gas, for the elec- 
tron/positron gas the exact expressions are integrated 
numerically (i.e. no assumptions about the degree of de- 
generacy or relativity are made) and the result is stored 
in a table. A sophisticated, biquintic Hcrmite polyno- 
mial interpolation is used to enforce the thermodynamic 
consistency at interpolated values(33). The photon con- 
tribution is treated as blackbody radiation. The EOS 
covers the density range form 10 -10 < pY e < 10 11 g 
cm~ 3 (Y" e being the electron fraction 5 ) and temperatures 
from 10 4 to 10 11 K. 

2.3. Gravity 

The self-gravity of the fluid is calculated via a par- 
allel version of the binary tree described in (jj). The 
same tree is used to search for the neighbor particles 
that are required for the density estimate, Eq. (JT|), and 



the gradients in Eqs. 



. The gas acceler- 
holc is treated in 



and 

ation due to a (SchwarzschTld) blac. 
the Paczyhski-Wiita approximation (|20fl . This approach 
has been shown (JTJ) to yield accurate results for the ac- 
cretion onto non-rotating black holes. To avoid numeri- 
cal problems due to the singularity at the Schwarzschild 
radius the pseudo potential is smoothly extended in a 
non-singular way down to the hole(|24[l with an absorb- 
ing boundary placed at a distance of 3GMbh/c 2 from the 
black hole with mass Mbh- 



In the presence of electron-positron pairs it is given by 
where n _/n + are the number densities of elec- 



trons/positrons. 
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2.4. A minimal nuclear reaction network 

To address whether tidal compression can trigger a 
thermonuclear explosion, we need to evolve the nuclear 
composition and to include the feedback onto the gas 
from the energy released by nuclear burning. Run- 
ning a full nuclear network with hundreds of species for 
each SPH particle would be computationally prohibitive, 
therefore a "minimal" network designed to provide the 
accurate energy generation (fTUl |3"2~)| is used. It cou- 
ples a conventional a-network stretching from He to Si 
with a quasi-equilibrium-reduced a-network. The QSE- 
reduced network neglects reactions within small equilib- 
rium groups that form at temperatures above 3.5 GK to 
reduce the number of abundance variables needed. Al- 
though a set of only seven nuclear species is used, this 
network reproduces all burning stages from He-burning 
to nuclear statistical equilibrium accurately. For more 
details and tests we refer to (fTOjl . 

In the presence of nuclear reactions the energy produced 
(or consumed) by nuclear reactions is given by 



£nuc,a — Na -B-j 
3 



dt 



(9) 



where is Avogadro's number, Bj is the nuclear bind- 
ing energy of the nucleus j, Yj a — n^ a j its abun- 
dance and nj ta is the number density of species j. Again, 
the subscript a indicates that these quantities are evalu- 
ated at the position of particle a. 

Since the nuclear reaction and the hydrodynamic time 
scales can differ by many orders of magnitude, the net- 
work is coupled to the hydrodynamics in an operator 
splitting fashion. In a first step, Eqs. ^ , Q , ^ , ^\ are 
integrated forward in time via a MacCormack predictor- 
corrector scheme (|13p with individual time steps (j2U| ) to 
obtain new quantities at time t n+1 . In this step we ig- 
nore the nuclear source term in Eq. pL the result is 
denoted by ti™ +1 . This value has to be corrected for the 
energy release that occurred from f" to t n+1 : 
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where p a (t) « p a (t n ) + F ^S^{p Q (i" +1 ) - p a {t n )} and 
T a (t) w T a (t n ) has been used to integrate the abundances 
Yj ja via the implicit backward Eulcr method (the network 
integration is described in detail in(fTtJ]0. The final value 
for the specific energy at time t n+1 is given by 

< +1 =< +1 +e Q , n ^ n+1 . (12) 

Now the EOS is called again to make all thermodynamic 
quantities consistent with this new value Once the 

derivatives have been updated, the procedure can be re- 
peated for the next time step. 

For the hydrodynamic time step we use the minimum of 
several criteria. We use a force criterion and a combina- 
tion of Courant-type and viscosity-based criterion(16) 



and 
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where f a is the acceleration, w SjQ is the sound velocity 
and [i a b a quantity used in the artificial viscosity tensor 
(for its explicit form see (|T5|l ). To ensure a close coupling 
between hydrodynamics and nuclear reactions in regions 
where burning is expected, we apply two additional time 
step criteria. One triggers on matter compression, the 
other on the distance to the black hole: 



and 



Aicomp,a = -0.03/(V • V) a 



Ai bh , a = 0.03/^/GM bh /r3 ha . 
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The "desired" hydro time step of each particle is then 
chosen as 

Aides, a = 0.2 min(Ai/ ia , Atc.a, At 

comp,ai 

A« bh , a ). (17) 

How these desired time steps are transformed into the 
individual block time steps is explained in detail in(29). 
As a test, low-resolution versions of the production runs 
were run once with this time step prescription and once 
with halved time steps. No noticeable differences could 
be found. 

One may wonder whether neutrino emission could 
subduct substantial amounts of energy during the burn- 
ing phase. To test for this, we have plotted in Fig. [2] 
the ratio of the nuclear energy produced and the energy 
lost to neutrinos within a typical pericentre passage of 
1 s duration. For the neutrino emission pair annihila- 
tion, plasma, photoneutrino, bremsstrahlung and recom- 
bination processes were considered according to the fit 
formulae of (jTTj) as coded by F. Timmes 6 . For the condi- 
tions relevant to this study neutrino emission was never 
relevant. 

3. RESULTS 

As initial condition we set up a white dwarf on a 
parabolic orbit around the black hole with an initial 

separation of several tidal radii i? t id = -Rwd ( Mwd ) 
In a large set of simulations we explored black hole 
masses of 100, 500, 1000, 5000 and 10 000 M Q and 
white dwarf masses of 0.2, 0.6 and 1.2 M . Each time 
several "penetration factors", P — i?tid/-R P erb where 
i?peri is the pericentre separation, were explored. To be 
conservative, we set the initial white dwarfs to very low 
temperatures (To = 5 • 10 4 K). The numerical resolution 
varied from 500 000 to more than 4 • 10 6 SPH particles. 
In all cases we found explosions (nuclear energy release 
larger than the white dwarf gravitational binding en- 
ergy) whenever the penetration factors exceeded values 
of about 3. Further details of these simulations will be 
discussed elsewhere (28). 

For conciseness we focus here on one exemplary sim- 
ulation of a 0.2 M , pure He white dwarf (modeled 
with more than 4 • 10 6 SPH particles) and a 1000 
Mq black hole, see Fig. [3] The first snapshot (0.12 
minutes after the simulation start) shows the stage 
of maximum white dwarf compression at pericentre 
passage, in which the white dwarf is also severely 
compressed perpendicular to the orbital plane. The 
peak compression occurs at a spatially fixed point (seen 
as the density peak in Fig. [3] left). The white dwarf fuel 

6 http://cococubed.asu.edu 
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is fed with free-fall velocity vg = (2GMbh/i? P eri) 1 ^ 2 = 

1.6 • 10 5 kms _1 (A/ bh /lOOOM ) 1/2 (i? pcri /10 4 km) _1/2 
into this compression point. The comparison with 
typical flame propagation speeds (~ 100 km/s) shows 
that combustion effects can be safely neglected for this 
investigation. 

During the short compression time (of order one second) , 
the peak density increases by more than one order of 
magnitude (with respect to the initial, unperturbed 
star) to > 6 • 10 5 g cm" 3 , the peak temperatures get 
close to nuclear statistical equilibrium (> 3.9 • 10 9 K). 
During this stage 0.11 M Q are burnt, mainly into silicon 
group (74%), iron group (22.5%) elements and carbon. 
The nuclear energy release triggers a thermonuclear 
explosion of the white dwarf. Since a much smaller 
fraction ends up in iron group nuclei (nickel), this 
explosion is undcrluminous in comparison to a normal 
type I a supernova which produces ~ 0.5 M Q of nickel. 
Due to the very different geometry, the lightcurves of 
such explosions are expected to deviate substantially 
from standard type la light curves. This topic deserves 
further detailed investigations. 

About 35 % of the initial stellar mass remain gravita- 
tionally bound to the black hole and will subsequently 
by accreted. During infall, matter trajectories become 
radially focused towards the pericentre. The large spread 
in the specific energy across the accretion stream width 
produces a large spread of apocentric distances and 
thus a fan-like spraying of the white dwarf debris after 
pericentre passage. This material interacts with the in- 



falling material in an angular momentum redistribution 
shock, see Fig. |4j which results in the circularization of 
the forming accretion disk. The subsequent accretion 
onto the black hole produces a soft X-ray flare close the 
Eddington-luminosity for a duration of months. 
For other combinations of black hole and white dwarf 
masses with P > 3 we found a similar behavior. 
Therefore, an underlumious thermonuclear explosion 
accompanied by soft X-ray flare may whistle-blow the 
existence of intermediate mass black holes in globular 
clusters. The rate of this particular type of explosion 
amounts to a few tenths of a percent of "standard" type 
la supernovae. Future supernova surveys such as SNF 
could be able to detect several of these events. 
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Fig. 1. — Sod shock tube problem in ID: the exact solution is given by the solid line, the numerical result is shown by the circles. In 
the lower left panel (pressure) the value of the time dependent artificial viscosity parameter a is overlaid. It is to be compared with the 
commonly used value a = 1 = const. 
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Fig. 2. — Comparison of nuclear energy produced (from an initial pure He composition) and the energy lost via various neutrino reactions 
during 1 s (the typical pericentre passage time, see Fig.pl. The thick solid line that starts at about log(T) = 9.1 and ends at about 
log(T) = 8.3 shows the trajectory of the hottest 10% of fne particles of the simulation shown in Fig. [3). To see clearly where neutrino 
emission becomes dominant, we have also thickened the contour where nuclear and neutrino contributions are equal. 
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Fig. 3. — Tidal disruption of a 0.2 Mq He white dwarf by a 1000 Mq black hole (located at the origin). 
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Fig. 4. — Zoom-in on the forming accretion disk. Color-coded is the column density (in g/cm 2 ), angular momentum is redistributed 
the shock that forms when the accretion stream interacts with itself. 



